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Abstract 

Given an approximation to a multiple isolated solution of a polyno- 
mial system of equations, we have provided a symbolic-numeric deflation 
algorithm to restore the quadratic convergence of Newton's method. Us- 
ing first-order derivatives of the polynomials in the system, our method 
creates an augmented system of equations which has the multiple isolated 
solution of the original system as a regular root. 

In this paper we consider two approaches to computing the "multiplic- 
ity structure" at a singular isolated solution. An idea coming from one of 
them gives rise to our new higher-order deflation method. Using higher- 
order partial derivatives of the original polynomials, the new algorithm 
reduces the multiplicity faster than our first method for systems which 
require several first-order deflation steps. 

We also present an algorithm to predict the order of the deflation. 
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1 Introduction 



This paper describes a numerical treatment of singular solutions of polynomial 
systems. A trivial example to consider would be a single equation with a double 
root, /(x) = = 0, or a cluster of two very close roots, f{x) = x'^ — = 0, 
where < £ ^ machine precision. In both cases getting good approximate solu- 
tions with straightforward numerical approaches such as Newton's method is not 
easy. Instead of attempting to solve the given equations we replace them with 
the system augmented by the equation's derivative, f{x) — {f{x),f'{x)) = 0. 
Note that this completely symbolic procedure leads to a system with exact reg- 
ular root in the first case, whereas in the second the system f{x) = would be 
inconsistent. However, a numerical solver applied to the latter would converge 
to a regular solution of a close-by system. 

In general setting, given an overdetermined system of equations in many 
variables with a multiple isolated solution (a cluster of solutions) our approach 
deflates the multiplicity of the solution (cluster) by applying a certain numerical 
procedure. From the point of view of the numerical analysis it may be called 
a reconditioning method: to recondition means to reformulate a problem so its 
condition number improves. 

Our deflation method was first presented at [25], and then described in 
greater detail in [T3] . In [H] , a directed acyclic graph of Jacobian matrices was 
introduced for an efficient implementation. 

On input we consider clusters of approximate zeroes of systems F(x) = 
{fi{x), f2{x), . . . , f]\[{x)) — oi N equations in n unknowns x G C". We 
assume the cluster approximates an isolated solution x* of F{x) — 0. Therefore, 
N > n. As X* is a singular solution, the Jacobian matrix of Fix), denoted 
by A{x), is singular at x* . In particular, we have r = Rank(A(cc*)) < n. 

In case r = n — 1, consider a nonzero vector A in the kernel of A{x*), which 
we denote by A G ker(A(a;*)), then the equations 

5^(-)=E^^^' ^ = 1,2,.. .,A^, (1) 

4=1 ^ 

vanish at a;*, because r = Rank(A(a;*)) < n. For r < n — 1, our algorithm 
reduces to the corank-1 case, replacing A{x) by A(x)B, where i? is a random 
complex 7V-by-(r -I- 1) matrix. For the uniqueness A G ker(A(a;*)), we add a 
linear scaling equation {h,X) = 1 (using a random complex (r -I- l)-vector h). 
and consider the augmented system 

( F{x) = 
G{x,X) = I A{x)BX = (2) 

I {KX) = 1. 

Let us denote by ^f{x*) the multiphcity of solution of the system 

F{x) — 0. In [2] we proved that there is a A* such that fia{x* , A*) < IjLp{x*). 
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Therefore, our deflation algorithm takes at most to — 1 stages to determine x* 
as a regular root of an augmented polynomial system. 

Related work. The literature on Newton's method is vast. As stated in [5], 
Lyapunov-Schmidt reduction (see also |6j §6.2], [T], [TT], and [H])) stands at 
the beginning of every mathematical treatment of singularities. We found the 
inspiration to develop a symbolic- numeric deflation algorithm in [19| . The sym- 
bolic deflation procedure of [T2] restores the quadratic convergence of Newton's 
method with a complexity proportional to the square of the multiplicity of the 
root. Algorithms to compute the multiplicity are presented in [5], [3], and [53]. 

Our Contributions. We establish the link between two different objects de- 
scribing what we call the multiplicity structure of an isolated singular solution: 
the dual space of differential functionals and the initial ideal with respect to a 
local monomial order, both associated to the ideal generated by the polynomials 
in the system in the polynomial ring. 

Next, following the latter method, we explain how to compute a basis of 
the dual space, flrst, following the ideas of Dayton and Zeng [3 , then using 
the approach of Stetter and Thallinger [23]. We provide a formal symbolic 
algorithm for each approach, respectively called the DZ and ST algorithms; the 
ingredients of the algorithms do not go beyond linear algebra. Moreover, we 
present an algorithm to determine the order of the deflation. 

The formalism developed for DZ and ST algorithms found a natural contin- 
uation in higher-order deflation method that generalizes and extends the first- 
order deflation in [14j . For the systems that require more than one deflation 
step by our first algorithm, the new deflation algorithm is capable of completing 
the deflation in fewer steps. 

Acknowledgements. The material in this paper was presented by the flrst 
two authors at the workshops on computational algebraic geometry and real- 
number complexity, organized respectively by Teresa Krick & Andrei Gabrielov 
and Peter Buergisser & Gregorio Malajovich. The authors thank the organizers 
for the opportunities to present their work at these FoCM 2005 workshops. 

2 Statement of the Main Theorem &; Algorithms 

The matrices A^'^'^ [x) we introduce below coincide for d = 1 with the Jacobian 
matrix of a polynomial system. They are generalizations of the Jacobian matrix, 
built along the same construction as the matrices used in the computation of 
the multiplicity by Dayton and Zeng in [3]. 

Definition 2.1 The deflation matrix ^('^^(a;) of a polynomial system F — (/i, 
/2; ■ • ■ , /w) of N equations in n unknowns x — (xi, X2, ■ ■ ■ , Xn) is a matrix with 
elements in C[x]. The rows of A'^'^^x) are indexed by x^fj, where \a\ < d 
and j — 1, 2, . . . , TV. The columns are indexed by partial differential operators 
= ^,f j where /3 7^ and < d. The element at row x"^ fj and 
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column of A'^'^^ (x) is 



■ 



(3) 



A^'^){x) has Nr rows and JVc columns, = N ■ ("+^"^) and A^c = - 1- 

Example 2.2 (Second-order deflation matrix) Consider a system of 3 equa- 
tions in 2 variables F = (/i,/2,/3) = 0, where fi = xf, f2 = — X2, and 
/a = x|. Then the second-order deflation matrix A^^\x) of / is 





dxt 


9x2 


92 


B 8 

^Xi ^X2 


52 

^X2 


/l 


- 2x1 





2 








/2 


2x1 


— 3^2 


2 





—6x2 


/3 





4xi 








12x1 


■Tl/l 


3x2 





6x1 








Xlf2 


3x2 




6x1 


3x2 


—6x1x2 


xifa 


x\ 


Ax 1X2 











X2fl 


2X\X2 


xl 


2X2 


2x1 





X2h 


2X\X2 


-4xi 


2X2 


2x1 


-12x^ 


X2h 













20xi 



(4) 



Notice that A^^^^ (x) (or the Jacobian matrix of F) is contained in the first three 
rows and two columns of A^'^\x). 

Definition 2.3 Let x* be an isolated singular solution of the system F{x) = 
and let d be the order of the deflation. Take a nonzero A^,.-vector {Xf})i3^o. \f3\<d 
in the kernel of ^(''^(a;*). It corresponds to what we call a deflation operator - 
a linear differential operator with constant coefficients A/3 



Q = 



E 

07^0, \0\<d 



We use Q to define Nr new equations 



(a:"/,) = 0, J = 1,2, 



,iV, \a\ < d. 



(5) 



(6) 



When we consider A/3 as indeterminate, we write gj^a as gj,a{x, A). In that case, 
we define m additional linear equations, for m = corank (^('^^(a;*)): 



'ifc(A) = ^6fc,/3A/3 - 1 = 0, k = l,2,...,m, 

where the coefficients bk,0 are randomly chosen complex numbers. 
Now we are ready to state our main theorem. 



(7) 
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Theorem 2.4 Let x* £ C" be an isolated solution of F{x) = 0. Consider the 
following system in C[x, A]; 

r f,{x) = 0, J = 1,2,..., TV; 
G^^Hx^X)^! g,Ax,\) = 0, j = l,2,...,7V, |a| <d; (8) 
[ hk{\) = 0, A: = l,2,...,m. 

For a generic choice of coefficients bk.p, there exists a unique A* e such 
that the system G^'^\x,\) has an isolated solution at (a;*. A*). Moreover, the 
multiplicity of (a;*. A*) in G{x, A) = is strictly less than the multiplicity of x* 
in F{x) = 0. 

To determine the order d, we propose Algorithm 12.51 This d is then used in 
Algorithmic 

Algorithm 2.5 d = MinOrderForCorankDrop(F. Xq) 

Input: _F is a finite set of polynomials; 

x'^ K X* , X* is an isolated multiple solution of F{x) = 0. 

Output: d is the minimal number such that the augmented system G*^''^ 
produced via a generic deflation operator Q of order d has 
corank of the Jacobian at x* lower than corank^(a;*). 

take a generic vector 7 = (71, . . . ,7„) e ker^(a;"); 
let H{t) = F{x° + jt) = F{x\ + -fit, ...,xl+ 7„t); 
d := min{ a \ a belongs to the support of H{t) } — 1. 

Algorithm 2.6 D'-'^^F = Deflate {F,d,XQ) 

Input: is a finite set of polynomials in C[a;]; 

d is the order of deflation; 

x'^ K X*, X* is an isolated multiple solution of F(x) = 0. 
Output: D'^'^^F is a finite set of polynomials in C[x, A] 

such that there is A* with ^]j{d)p{x* , \*) < ^p{x*). 

determine the numerical corank m of A'^'^\x'^); 
return D^F ■= G^'^\x,\) as in ©. 

3 Multiplicity Structure 

This section relates two different ways to obtain the multiplicity of an isolated 
solution, constructing its multiplicity structure. Note that by a "multiplicity 
structure" - a term without a precise mathematical definition - we mean any 
structure which provides more local information about the singular solution in 
addition to its multiplicity. In this section we mention two different approaches 
to describe this so-called multiplicity structure. 
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Example 3.1 (Running example 1) Consider the system 



F{x) 



x\x\ = 

xf + x\x2 — 0. 



(9) 



The system F{x) — has only one isolated solution at (0, 0) of high multiplicity. 



3.1 Standard Bases 

Assume G C" is an isolated solution of the system F{x) = 0. Let / = (F) C 
R = C[x] be the ideal generated by the polynomials in the system. Given a local 
monomial order >, the initial ideal in>(J) — {in>(/) \ f £ 1} C R describes the 
multiplicity structure of by means of standard monomials, i.e.: monomials 
that are not contained in in> (/) . A graphical representation of a monomial 
ideal is a monomial staircase. 

Example 3.2 (Initial ideals with respect to a local order) Consider the 
system © of Example 13.11 

Figure [1] shows the staircases for initial ideals of / = (F) w.r.t. two local 
weight orders >i^. Computer algebra systems Macaulay 2 [7] and Singular [TU] 
can be used for these kind of computations, see also [ll[9j for theory, in partic- 
ular, on Mora's tangent cone algorithm [T7j. 

In the example the leading monomials at the corners of the staircase come 
from the elements of the corresponding standard basis. For the weight vector 
w = (—1, —2) the original generators give such a basis (initial terms underlined). 
For w — (—2,-1) one more polynomial is needed. o 



Figure 1: Two monomial staircases for two different monomial orderings applied 
to the same system. The full circles represent the generators of the initial ideals. 
The multiplicity is the number of standard monomials, represented by the empty 
circles under the staircase. 



Below we will show how to compute the multiplicity of (0, 0). 



o 
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3.2 Dual Space of Differential Functionals 

Another approach at the multiphcity structure is described in detail in [211 [H] ; 
see also [H] . Using duality to define the multiplicity of a solution goes back to 
Macaulay jT6j . In this approach, differential Junctionals are denoted by 



1 (9l"lf 
A«(/) ^ 



ai! •■•«„! dx"^---dx°'" 
Observe that 



(10) 



x=0 



We then define the local dual space of differential functionals Dq [I] as 

Do[I] = {Le Span{A„ | a e Z|o}| L{f) = for all / e /}, (12) 

Example 3.3 (Dual space of running example 1) For the ideal defined by 
the polynomials in the system ([5]) we have 

I?o[/] = Span{ A(4^o) - A(3^i), A(3,o) , ^(2,1) , ^(1,2) , 

A(2,o) , 1) , A(o 2) , A(i 0) , A(o 1) , A(n,o) }• 

Notice that here the basis of the dual space is chosen in such a way that the 
(underlined) leading terms with respect to the weight order >(2,i) correspond 
to the monomials under the staircase in Example l3.1l for the order >(-2,-i)- We 
will show that it is not a coincidence later in this section. o 



3.3 Dual Bases versus Standard Bases 

Since both local dual bases and initial ideals w.r.t. local orders describe the 
same, there exists a natural correspondence between the two. 

Let > be an order on the nonnegative integer lattice Z"q that defines a local 
monomial order and let ^ be the opposite of >: i.e. a ^ P ^ a < f3. (Note: cz 
defines a global monomial order.) 

For a linear differential functional L = ^ Cq A^ define the support: supp(L) = 
{a e Z>Q I Ca ^ 0}. For the dual space, supp(-Do[/]) — UlgUoI/] supp(-^)- 

Using the order >z we can talk about the leading or initial term of L: let 
in^(i) be the maximal element of supp(L) with respect to Define the initial 
support of the dual space as my{Do[I]) — {in^(i) | L G Do[I]}. The ini- 
tial support is obviously contained in the support, in our running example the 
containment is proper: 

in(2,i)(i^o[/]) = <3}U{{A,0)} 

C {(ij) M+j<3}U{(4,0)U(3,l)} = supp(i?o[/]). 
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Theorem 3.4 The number of elements in the initial support equals the dimen- 
sion of the dual space, therefore, is the multiplicity. Moreover, with the above 
assumptions on the orders > and y, the standard monomials w.r.t. the local 
order > are {x" \ a G in^(_Do[-^])}- 

Proof. Pick e Dq[I],/3 E my{DQ[I]) such that my{Lp) = (3. One can easily 
show that {Lp} is a basis of -Do[-^]- 

Take a monomial G in>(/), then there is / G / such that x" — in>(/). 
Next, take any linear differential functional L with in^(L) = a. Since the orders 
> and y are opposite, there are no similar terms in the tail of L and the tail of 
/, therefore, L(/) = in^(L)(in>(/)) ^ 0. 

It follows that, L ^ -Do[^], which proves that the set of standard monomials 
is contained in the initial support of -Do[^]- They are equal since they both 
determine the dimension. □ 

Consider the ring of linear differential operators P = C[c?] with the natural 
action (denoted by "■") on polynomial ring R = C[x]. 

Lemma 3.5 Let Q G C[c?] and f G C[a;] such that in^(Q) >: in>(/) (in Z"gJ. 
Then in>(Q • /) = in>(/) - in^(O) G Z^^- 



4 Computing the Multiplicity Structure 

Let the ideal / be generated by /i, /2, . . . , /at. Let D^q \i] the part of Do[I] 
containing functionals of order at most d. We would like to have a criterion 
that for the differential functional L of degree at most d guarantees L G [I] . 

Below we describe two such criteria referred to as closedness conditions; 
their names are arranged to match the corresponding computational techniques 
of Dayton-Zeng [3^ and Stetter-Thallinger [23] that we will describe later re- 
spectively as the DZ and ST algorithms. 

A functional i = Ca Aq, with Cq, G C of order d belongs to the dual space 
Do[I] if and only if 

• (DZ-closedness) L{g ■ fi) = for all i = 1, 2, . . . , and polynomials 
g{x) of degree at most d — 1. 

• (ST-closedness) L{fi) — for all i and o'j{L) G Dq[I] for all j = 
1, 2, . . . , n, where ctj : 13o[-^] ~^ ^o[^] is a linear map such that 

-.(A.) = (.'^' (14) 
■' y ^a-cj , otherwise. ^ ' 

The basic idea of both DZ and ST algorithms is the same: build up a basis 
of Do incrementally by computing D^q ^ for d ^ 1,2,... using the corresponding 
closedness condition. The computation stops when D^q ^ ~ D^q ^\ 
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Example 4.1 (Running example 2) Consider the system in C[xi,a;2] given 
by three polynomials /i — xiX2, f2 — xl — X2, and /s — which has only one 
isolated root at (0,0). o 



4.1 The Dayton-Zeng Algorithm 

We shall outline only a summary of this approach, see [3] for details. 
If is a solution of the system, then D^^'' — Span{Ao}. 
At step 0? > 0, we compute Let the functional 

L= J2 c„A„ (15) 

belong to the dual space . Then the vector of coefficients Ca is in the kernel 
of the following matrix m|^^ with NB{d — 1) rows and B{d) — 1 columns, where 
i3((i) — ^^f) is the number of monomials in n variables of degree at most d. 

The rows of iV/jj'^ are labelled with a;"/^, where jaj < d and j = 1, 2, . . . , TV. 
The columns correspond to A^j, where /3 ^ 0, < d. 

[The entry of Afj^''^ in row 2^"/^ and column A^] — Ap{x" fj). (16) 



At the step d = 3 we have the following AI- 





A(i,o) 


A(o,i) 


A(2,0) 


^(1.1) 


^(0,2) 


^(3,0) 


A(2,l) 


A(l,2) 


A{0,3) 


fl 











1 

















/2 








1 





-1 














/3 





























Xlfl 




















1 








Xlf2 

















1 





-1 





Xlfs 





























X2fl 























1 





2:2/2 




















1 





-1 


2:2/3 





























x'ifi 






























Note that the last block of 9 rows is entirely zero. 

Analyzing the kernel of this matrix one sees that there are no functionals of 

(2) 

degree 3 in the dual space, which is then is equal to Dq [I] 

Do[/] =Span{A(o,o),A(i,o),A(o,i)'^(2,o)+A(o,2)}. (17) 



4.2 The Stetter-ThaUinger Algorithm 

The matrix Mg^r^ is a matrix consisting of n + 1 blocks stacked on top of each 
other: 
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• The top block contains the first N rows of M^f^; 

Ad) 



For every j = 1, 2, . . . , n, let 5"] be the {B{d-1) - 1) x {B{d) - l)-matrix 
for the linear map 

a, : ni^^/ Span{Ao} ^ i?^,''"''/ Span{Ao} (18) 



w.r.t. standard bases of functionals. 

The block Mg"^ ^''Sj represents the closedness condition for the "anti- 
derivation" a-j. 



Let us go through the steps of the algorithm for the Example 14.1 
Step 1. At the beginning we have equal to 





^(1,0) 


^(0,1) 


/l 
















h 









Therefore, D^^^ = Span{A(o,o), A(i^o), A(o,i)}. 
Step 2. Since M^g^sf^ = for all j, the matrix M 



(2) 
ST 



IS 





A(i,o) 


A(o,i) 


A(2,0) 




A(0,2) 


/l 











1 





/2 








1 





-1 


h 



































)}■ 



Therefore, D}^ = Span{A(o,o) , A(i^o), A(o,i), A(2,o) + A(o,2 

(21 

We can "prune" the matrix Mgj, by row-reducing it to the following matrix 
with the same kernel: 



M 



(2) 
ST 








1 

-1 



(3) 

Step 3. Compute S\ that represents ai 





A(i,o) 


A(o,i) 


A(2,0) 


A(i,i) 


A{0,2) 


A(3,0) 


A(2,l) 


A{1,2) 


A(0,3) 


A(i,o) 








1 




















A(o,i) 











1 

















A(2,0) 

















1 











A(i,i) 




















1 








A(0,2) 























1 
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The matrix 5*2 can be defined similarly. 
The top block of the matrix Mg^ is 







A(0,1) 


^(2,0) 


^(1.1) 


^(0,2) 


^(3,0) 


A(2,l) 


A(l,2) 


A(0,3) 


/l 











1 

















h 








1 





-1 














h 






























Despite the last 4 columns being 0, there are no new elements of order 3 in 

~ (2) (3) 

the dual space due to the other two blocks: MgrfS{ : 





A(i,o) 


A(o,i) 


A(2,0) 


A(i.i) 


A(0,2) 


A{3,0) 


A(2,l) 


A{1,2) A(o,3) 


Xlfl 




















1 





Xlf2 

















1 





-1 



and M^^2^4' 





A(i,o) 


A(o,i) 


A(2,0) 


A(i,i) 


A(0,2) 


A{3,0) 


A(2,l) 


A{1,2) 


A{0,3) 

























1 





X2f2 




















1 





-1 



Comparing to DZ algorithm, in step 3, we managed to avoid the computation 

(3) 

of 9 last zero rows of Mj^^ in this particular example. We now also see how its 
4 last nonzero rows show up in the "closedness condition" blocks of Mjj^. 



5 Proofs and Algorithmic Details 

In this section we justify the main theorems stated before and give details about 
the algorithms presented above. 

5.1 First-Order Deflation 

In this section we summarize our deflation method introduced in |14| . Not only 
it is done for the convenience of the reader, but also for our own convenience 
as we plan to build a higher-order deflation algorithm in Section O using the 
algorithm following the pattern established in this section. 

One deflation step vifith fixed A. The basic idea of the method is relatively 
simple. Let A G C" be a nonzero vector in ker(^(a;*)), then the equations 

g^{x) = \-Vf,{x) = Y,X,^^, z = l,2,...,7V (19) 

OX j 

z— 1 

have X as a solution. Moreover, 
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Theorem 5.1 The augmented system 

G{x)^{h,...jN,9i....,9N){x)^0 (20) 

of equations in C[x\ is a deflation o/ the original system F{x) = at x* , i.e. 
G{x*) ~ and the multiplicity oj the solution x* is lower in the new system. 

The original proof of this statement in .H] uses the notion of a standard 
basis of the ideal / = (/i, /2, ■ • ■ , /jv) in the polynomial ring R = C[x] w.r.t. 
a local order; this tool of computational commutative algebra can be used to 
obtain the multiplicity of x*, which is defined as the C-dimension of the local 
quotient ring Rx/Rxl- 

On the other hand it is in correspondence with another way of looking at 
multiplicities - dual spaces of local functionals, so the proof can be written in 
that language as well (see Section [S]) . 

One deflation step with indeterminate A. Without loss of generality, 
we may assume corank (A(a;*)) — 1; consult fT^ to see how the general case 
is reduced to this. Consider + 1 additional polynomials in C[a;,A] in 2n 
variables: 

g,{x,\)^\-VMx)^J2^^^^' (* = 1,2,...,7V) (21) 

j=l 

n 

h{\)=Y.^,X,-l, (22) 

where the coefficients bj are random complex numbers. 

Theorem 5.2 Let x* G C" be an isolated solution of F(x) = (in C[x]). 

For a generic choice of coefficients bj , j = 1,2, ... , n, there exists a unique 
A* G C" such that the system 

G{x, X)^{f,,...jN,gi,...,gN,h){x,\) = (23) 

of equations in C[x,\] has an isolated solution at (a;*, A*). 

The multiplicity of (a;*, A*) in G{x,X) ~ is lower than that of x* in 
F{x) = 0. 

Proof. Follows from Proposition 3.4 in [T3]. □ 

Theorem 15.21 provides a recipe for the deflation algorithm: one simply needs 
to keep deflating until the solution of the augmented system corresponding to 
X* becomes regular. 

As a corollary we have that the number of deflations needed to make a 
singular isolated solution x* regular is less than the multiplicity of x* . 
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5.2 Higher-Order Deflation with Fixed Multiphers 

We use the deflation operator to define an augmented system. 

Theorem 5.3 Let /i, /2, . . . , /at form a standard basis of I w.r.t. the order 
opposite to ^. Consider the system G^'^'^{x) = in C[x], where 

I 9] Ax) [j ^1,2,...,N, \a\<d) ^ ^ 

(a) The system G^'^\x) ~ is a deflation of the original system F{x) = 
at X*. 

(b) Let I = (F) and J = (G^'*-') be the ideals generated by polynomials of the 

systems and > be a global monomial order on Z>q. Then the following 
relation holds for initial supports 

in^{Do[J]) ^ {P - Pq I /3Gin^(i?o[/])}nZ|o, (25) 

where (3q is the maximal element of the set in'^{Do[L\) fl {(3 : \[3\ < d}. 

Proof. Let A e kcT{A'-'^'> {x*)) be the vector used above to construct the operator 
Q e C[d] and the equations gj_a{x) = 0. 

First of aU, gj^a{x*) — {Q ■ {x" fj))\x=x' = provided |a| < d (by construc- 
tion), hence, x* is a solution to G^'^^x) = 0. 

To prove (a), it remains to show that the multiplicity drops, which follows 
from part (b) that is treated in the rest of this proof. 

We shall assume for simplicity that x* — Q. This is done without the loss of 
generality using a linear change of coordinates: a; i— > a; + a;*. It is important to 
note that in the new coordinates polynomials Q ■ [x"- fj{x + x*)) generate the 
same ideal as the polynomials Q ■ {{x — x*)" fj{x + x*)). 

Recall that L = (F) = (/i, /2, . . . , /at), let J = (G^'')) D / be the ideal gen- 
erated by the polynomials in the augmented system. The reversed containment 
holds for the dual spaces: Do[I] D Dq[J]. 

There is a 1-to-l correspondence between linear differential operators and 
linear differential functionals: 

Let (f) : C[d] Dq and t : Dq — > C[d] be the corresponding bijections. 

As in Section[3]we order terms with ^, a global monomial order. Notice 
that since the choice of coefficients of the operator Q is generic, /3q = in^ (Q) = 
iiiy{(j>{Qj) is the maximal element of the set iny{Do[L]) n {(3 : |/3| < d}. 

Next, we use the condition that fi form a standard basis. Since the corners 
of the staircase correspond to the initial terms of fi, by Lemma 13.51 the staircase 
created with the corners at in>(Q • [x'^fi)) bounds the set {P — Pq \ fi € 
in^(Z?o[-f])} n Z"q, which, therefore, contains the initial support of Uoi"/]- D 
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Corollary 5.4 // there exist a local monomial order > such that the minimal 
(standard) monomial in the set {x" ^ in>(/) : \a\ < d} is minimal in the set of 
all standard monomials, then x* is a regular solution of G^'^\x) = 0. 

Ideally we would like to be able to drop the assumption of the original 
polynomials forming a standard basis, since computing such a basis is a complex 
symbolic task, whereas our interest lies in the further numericalization of the 
approach. The following weaker statement works around this restriction. 

Assuming x* = 0, let supp(F) = Uj=i,2,...,Ar supp(/j). 

Proposition 5.5 Assume A{0) = 0. Let do = min{|Q;| : G supp(F)}. 

Then, in the notation of Theorem \5.!A for a generic deflating operator Q 
the system G^'^\x) — 0, where, d < do is a deflation of the original system 
F{x) = at the origin. 

Moreover, if d = do — 1 then the Jacobian of G^''-' (0) is not equal to zero. 

Proof. Fix a local monomial ordering that respects the degree. With the above 
assumptions, the initial ideal in((_F)) will contain monomials of degree at least 
do- On the other hand, for a generic choice of the deflating operator Q the 
support supp(G('')) would contain a monomial of degree less than do- Therefore, 
there exists a monomial in in ((supp(G('')))) that is not in in ((F)), hence, G^'*) 
is a deflation. 

If d = |(io| — 1, then there is such monomial of degree 1, which means that 
the Jacobian of the augmented system is nonzero. □ 

Remark 5.6 Note that if the deflation order d is as in Proposition [531 then it 
suffices to take an arbitrary homogeneous deflation operator of order d. 

Next we explain the practical value of Proposition 15.51 Let K = kerA(O) 
and c = corankA(O) — dimi^. Without a loss of generality we may assume 
that K is the subspace of C" has {xi, . . . ,Xc} as coordinates. 

Now consider the system F'{xi, . . . ,Xc) = F{xi, . . . , Xc, 0, . . . , 0). This sys- 
tem has an isolated solution at the origin, and Proposition 15.51 is applicable, 
since the Jacobian is zero. Moreover, if we take the deflation of order d = do — l 
of the original system F, with do coming from the Proposition, the corank of 
the Jacobian the augmented system G^"^^ is guaranteed to be lower than that 
of A(0). 

Let us go back to the general setup: an arbitrary isolated solution x*, the 
Jacobian A{x*) with a proper kernel K, etc. Algorithm 12.51 is a practical algo- 
rithm that can be executed numerically knowing only an approximation to x*. 

Proof. [Proof of correctness of Algorithm 12.51 for a;° = x*.] We can get to 
the special setting of Proposition 15.51 in two steps. First, apply an afline trans- 
formation that takes x* to the origin and kei A{x*) to the subspace K of C" 
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spanned by the first c = corankA(£c*) standard basis vectors. Second, make a 
new system F'{xi, . . . , Xc) = by substituting the Xi = in F for i > c. 

Let 7' = (7^, . . . ,7^) e K he the image of the generic vector 7 under the 
hnear part of the affine transform. Then H{t) ~ F'{'~f[t, . . . , j'^t). 

Since 7' is generic, the lowest degree do of the monomial in supp(-F') is equal 
to min{a | g suppi?(i)}. According to the Proposition 15.51 and the discussion 
that followed, d = do — 1 is the minimal order of deflation that will reduce the 
rank of the system. □ 

Remark 5.7 In view of Remark l5.6l it would be enough to use any homogeneous 
deflation operator of order d: 



such that the vector A of its coefficients is in the kernel of the truncated deflation 
matrix, which contains only the rows corresponding to the original polynomials 
F and only the columns labelled with with |/3| = d. 

5.3 Indeterminate Multipliers 

As in Section l5.1i wc now consider indeterminate A/3. Now we should think of 
the differential operator L{\) e C[A, d] and of additional equations gj^aix, A) G 
C[x, A] as depending on A. 

Proof of Theorem \2.4\ Picking m — corank (A'^''' (a;*)) generic linear equations 
hk guarantees that for x = x* the solution for A exists and is unique; therefore, 
the first part of the statement is proved. 

The argument for the drop in the multiplicity is similar to that of the proof 
of Theorem O □ 

6 Computational Experiments 

We have implemented our new deflation methods in PHCpack [5S] and Maple. 
Below we report on two examples. 

One crucial decision in the deflation algorithm is the determination of the 
numerical rank, for which we may use SVD or QR in the rank-revealing algo- 
rithms. Both SVD and QR arc numerically stable. We summarize the result 
from [H page 118] for the problem of solving an ovcrdetermined linear sys- 
tem Ax — b. The solution obtained by QR or SVD minimizes the residual 
IK^l -I- SA)x — {b + Sb)\\2 where the relative errors have the same magnitude as 
the machine precision e: 




(27) 



|/3|=d 




(28) 
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6.1 Running Example 1 



To find initial approximations for the roots of the system we must first 
make the system "square", i.e.: having as many equations as unknowns, so we 
may apply the homotopies available in PHCpack ^5\. Using the embedding 
technique of [20] (see also [21]), we add one slack variable z to each equation of 
the system, multiplied by random complex constants 71, 72, and 73: 



Observe that the solutions of the original system F{x) = occur as solutions of 
the embedded system E{x, z) = with slack variable z ~ 0. At the end points 
of the solution paths defined by a homotopy to solve E{x,z) = 0, we find nine 
zeroes close to the origin. These nine approximate zeroes are the input to our 
deflation algorithm. 

The application of our first deflation algorithm in [Tl] requires two stages. 
The Jacobian matrix of F{x) — has rank zero at (0, 0). After the first defiation 
with one multiplier, the rank of the Jacobian matrix of the augmented system 
G'(a;,Ai) = equals one, so the second deflation step uses two multipliers. 
After the second deflation step, the Jacobian matrix has full rank, and (0, 0) 
has then become a regular solution. Newton's method on the final system then 
converges again quadratically and the solution can be approximated efficiently 
with great accuracy. Once the precise location of a multiple root is known, we are 
interested in its multiplicity. The algorithm of [3] reveals that the multiplicity 
of the isolated root equals seven. 

Starting at a root of low accuracy, at a distance of 10~^ from the exact 
root, the numerical implementation of Algorithm 12.51 predicts two as the order, 
using 10~^ as the tolerance for the vanishing of the coefficients in the univariate 
interpolating polynomial. The Jacobian matrix of the augmented system G*-^^ 
has full rank so that a couple of iterations suffice to compute the root very 
accurately. 

6.2 A Larger Example 

The following system is copied from [12] : 



Counted with multiplicities, the system has 54 isolated solutions. We focus on 
the solution (0,0, —1) which occurs with multiplicity 18. 

Although Algorithm 1 suggests that the first-order deflation would already 
lower the corank of the system, we would like to search for a homogeneous 
deflation operator Q of order two. 




(29) 




2xi + 2x1 + 2x2 + 2x^ + a;§ - 1 = 
{xi + a;2 - - 1)^ - = 
{2x1 + 2a:l + 10x3 + bxj + 5)^ - lOOOx^ = 0. 



(30) 
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To this end we construct the (truncated) deflation matrix A{xi, X2, X3) with 
has 12 rows and only 6 columns, which correspond to {df, 8182,8183, 9|, 8283, 83}. 

The kernel of A(0, 0, -1) is spanned by (1, 6, 8, -3, 0, 4)'^ and (0, 3, 3, -1, 1, 2)'^. 
The operator corresponding to the former, 

Q = 8f + 68182 + 88183 - 39| + A8i, (31) 

regularizes the system, since the equations 

Q ■ {xifi) = 8x1 + 24x2 + 16x3 + 16 = 
g-(x2/i) - 24x1 - 24x2 = (32) 

Q-{x3h) = 32x1 + 16x3 + 16 = 0. 

augmented to the original equations, give a system with the full-rank Jacobian 
matrix at (0, 0, —1). 

7 Conclusion 

In this paper we have described two methods of computing the multiplicity 
structure at isolated solutions of polynomial systems. We have developed a 
higher-order deflation algorithm that reduces the multiplicity faster than the 
flrst-order deflation in [14j. 

In our opinion, one of the main benefits of the higher order deflation for 
the numerical algebraic geometry algorithms is the possibility to regularize the 
system in a single step. For that one has to determine the minimal order of such 
a deflation or, even better, construct a sparse ansatz for its deflation operator. 
Predicting these numerically could be a very challenging task, which should be 
explored in the future. 
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